Mid- J CO Emission From NGC 891: Microturbulent Molecular 
Shocks in Normal Star Forming Galaxies 
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ABSTRACT 

We have detected the CO(6-5), CO(7-6), and [C I] 370 ^m lines from the nuclear region of 
NGC 891 with our submillimeter grating spectrometer ZEUS on the CSO. These lines provide 
constraints on photodissociation region (PDR) and shock models that have been invoked to 
explain the H 2 S(0), S(l), and S(2) lines observed with Spitzer. We analyze our data together 
with the H2 lines, CO(3-2), and IR continuum from the literature using a combined PDR/shock 
model. We find that the mid-J CO originates almost entirely from shock-excited warm molecular 
gas; contributions from PDRs are negligible. Also, almost all the H2 S(2) and half of the S(l) 
line is predicted to emerge from shocks. Shocks with a pre-shock density of 2 x 10 4 cm~ 3 and 
velocities of 10 km/s and 20 km/s for C-shocks and J-shocks, respectively, provide the best fit. 
In contrast, the [C I] line emission arises exclusively from the PDR component, which is best 
parameterized by a density of 3.2 x 10 3 cm -3 and a FUV field of G = 100 for both PDR/shock- 
type combinations. Our mid-J CO observations show that turbulence is a very important heating 
source in molecular clouds, even in normal quiescent galaxies. The most likely energy sources for 
the shocks are supernovae or outflows from YSOs. The energetics of these shock sources favor 
C-shock excitation of the lines. 

Subject headings: turbulence — galaxies: ISM — galaxies: star formation — galaxies: individual: NGC 
891 — submillimeter: ISM — techniques: spectroscopy 



1. Introduction 

Turbulence and shocks play a profound role in 
the process of star formation and nuclear gas ac- 
cretion. Shocks compress the interstellar medium 
(ISM) and can trigger star formation, while turbu- 
lence prevent the collapse of the ISM. These mech- 
anisms can also remove angular momentum from 
the gas, leading to gas infall toward galactic nuclei, 
a process that is important in interacting galaxies. 

Both mechanisms are very important in con- 
verting mechanical energy into thermal energy and 
can dominate the heating budget of molecular 
clouds. The main coolants of the shock heated 
molecular gas are H2 and CO lines. With their 



improved sensitivity, recent space observatories 
(ISO, Spitzer) were able to trace extragalactic ro- 
tational H2 emission and revealed the importance 
of this emission and the prevalence of shock ex- 
cited molecular gas. For example, H2 emission 
associated with shock or microturbulence excited 
molecular gas was found in the group Stephan's 
Quintet ( |Appleton et al J2006| |Cluver et al.|2010[ ), 
the radio galaxy 3C 326 ( Ogle et al.|2007 ), nearby 
edge-on ga laxies NGC 4 565, NGC 5907 (jLaiile 
et al.|2010[ ) and NGC 891 ( |Stacey et al.|2010[ here- 



after SCB10), as well as in the extended cold neu- 
tral medium in the Milky Way (Falgarone et al. 



2005). Ground-base mid-J CO observations are 
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also good probes of shock or turbulence excited 
molecular gas. Bradford et al. (2005) have in- 
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voked dissipation of MHD turbulence to explain 
their mid-J CO observations of the circumnuclear 
ring in Galactic center. In the nuclear region 
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of NGC 253, t he strong mid- J 12 CO (e.g. [Brad- 
ford et al.|2003|) and 13 CO(6-5) (jHailey-Dunsheath 



et al.||2008 ) emission is consistent with heating of 
the warm molecular gas by either decay of super- 
sonic turbulence through shocks or by an enhanced 
cosmic ray density. 

Due to its proximity (9.5 Mpc), and edge-on 
presentation (i > 89°; |Sancisi fc Allen|1979| , NGC 
891 is the primary target for studies of the effects 
of energetic processes associated with star forma- 
tion on the interstellar medium in normal spiral 
galaxies. The edge-on view of NGC 891 results 
in greatly enhanced gas columns so that relatively 
weak lines such as the low-lying quadrupole rota- 
tional transitions of H2 become observable. The 
mid-IR H2 rotational transitions in NGC 891 were 
studied for the first time with ISO (iValentijn & 



van der Werf 1999) and then more recently with 



Spitzer (SCB10). The latter study combined the 
H2 lines observed with Spitzer (S(2), S(l), and 
S(0)) together with [C II], [O I] and IR contin- 
uum observations from ISO to constrain gas exci- 
tation mechanisms. They found that most of the 
S(0) line emission arises in photodissociation re- 
gions (PDRs), and most of the S(2) line emission 
arises in low velocity microturbulent shocks, with 
the S(l) line arising from both types of sources in 
the ratio 70% PDRs, and 30% shocks. 

Here, we expand on the analysis of |SCB10| by 
including submillimeter observations of the CO (6- 
5), CO(7-6), and [C I] 370 (jm lines. The mid-J 
CO lines can arise from gas heated by a variety of 
mechanisms, including FUV photons (PDR), X- 
rays (XDR), cosmic rays, and shocks. The [C I] 
emission arises mainly from PDRs. Since the 
Spitzer H2 study showed that the main heating 
sources of the warm molecular gas in NGC 891 
are turbulence and PDRs we will focus on these 
two heating mechanisms and use our new submil- 
limeter lines to both test and constrain the models 
of |SCB10[ As in their analysis, we apply the PDR 
model of Kaufman et al. ( 1999 ) and the shock 



model of Flower & Pineau Des Forets (20101 



2. Observations 

We have observed the 12 CO J = 6 — > 5 (433.56 
/tm) and J = 7 — > 6 (371.65 /xm) rotational tran- 
sitions and the [C I] 3 P 2 ^ 3 Pi (370.41 /on) 
fine structure line with our submillimeter grating 



spectrometer ZEUS ( Hailey-Dunsheath |2009 ) on 
the Caltech Submillimeter Observatory (CSO) on 
Mauna Kea, Hawaii, in December 2006. 

ZEUS employs a 1 x 32 thermistor sensed detec- 
tor array that is optimized for observations in the 
submillimeter wavelength regime. It provides 32 
spectral elements at one spatial position. For the 
Dec. 2006 observing run, we mounted two band- 
pass filter appropriate for the 350 and 450 /xm 
telluric window adjacent to each other at the en- 
trance slit of the detector housing, each filter cov- 
ering 16 spectral pixels. This allowed us to observe 
in the 350 and 450 /im telluric windows simulta- 
neously and it still provided a full bandwidth cov- 
erage of about 5000 km/s in each band. The spec- 
tral resolution of ZEUS is A/AA ~ 1000 across the 
350 /im and 450 /im telluric windows. The CO(7- 
6) and [C I] lines, which are separated in velocity 
by only 1001 km/s, were observed simultaneously 
within a single spectrum. We verified the spectral 
calibration through observation of CO gas absorp- 
tion lines from a calibration unit mounted in front 
of the ZEUS dewar. 

The CO lines were observed on two different 
nights. On the first night we observed only the 
CO(6-5) emission at the position RA=02 h 22 m 33.10 s 
and Dec=+42°20'55.1" (J2000). On the second 
night we observed the CO (6-5), (7-6), and [C I] 
lines, but centered ~ 2" further to the north. This 
difference is, however, within the pointing accu- 
racy of the telescope. All observations were made 
in chop-nod mode with the secondary chopping at 
2 Hz and a chop throw of 30". 

We established a pointing model from observa- 
tions of Saturn and Uranus at various elevation 
and azimuth positions. The derived pointing ac- 
curacy is < 5". We estimated the coupling of the 
planets to the beam by assuming the planets to 
be disk-like sources with uniform emissivity and 
assuming the ZEUS beam to be Gaussian. During 
the observations the beam size was 10" in the 350 
/*m band and 11" in the 450 fim band. At the 
distance of NGC 891, 10" corresponds to 460 pc. 

All observations were corrected for atmospheric 
absorption through use of the CSO radiometers. 
Each spectrum was corrected for the response 
function of the grating and detector array, then 
flux calibrated with respect to Saturn, which was 
assumed to emit like a blackbody with a tempera- 



ture of T = 113 K at 370 /im and 434 /im (Hilde- 
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brand et al. 1985). The correction for telluric 



transmission is the largest source of error. Overall, 
we estimate a line flux uncertainty of ~ 25%. 

The CO(6-5), CO(7-6), and [C I] 370 fim lines 
were all strongly detected (Fig. [I]) , and their ob- 
served line intensities are listed in Table [l] 

3. Origin Of The Mid-J CO Emission in 
NGC 891 

The warm molecular gas traced by mid-J CO 
emission can be excited by several mechanisms. 
For spiral galaxies with quiescent star formation, 
the most relevant are shocks in micro-turbulent 
gas or far-UV radiation in photodissociation re- 
gions (PDRs) . Spitzer observations of the H2 mid- 
IR rotational transitions in NGC 891 have revealed 
the presence of warm, shock-excited molecular gas. 
The line intensities can only be explained if part 
of the H 2 emission arises from this warm, micro- 



turbulent gas phase ((SCB10|. While the H 2 S(0) 
likely originates entirely from PDRs, 30% of the 
H 2 S(l) and the majority (80%) of the H 2 (S2) 
emission likely emerges from gas excited by C-type 
shocks with a pre-shock gas density of ~ 2 x 10 4 
cm -3 , shock velocities between 20-30 km/s, and a 



shock filling factor of 2.8 (SCBlOl. The properties 



of this shock excited gas should also be imprinted 
in the line ratio of our CO observations. 

We are taking a slightly different approach in 
our shock-analysis than |SCB10[ For their analysis 
they considered H 2 emission from within large ar- 
eas of the galaxy since they compared their Spitzer 
observations with far-IR lines obtained with ISO, 
which had a beam size of 75". Here, we focus 
on the central 11" of NGC 891 and with our new 
ZEUS mid-J CO observations we can refine the 
shock modeling. We combine the Spitzer H 2 ob- 
servations (Table |2| with our ZEUS observations 
and apply the shock models of |Flower fc Pineau] 
Des Forets ( 2010[ ). The spatial resolution of the 
Spitzer H 2 observations are very close to our ZEUS 
observations. The slit orientation of the long-high 
IRS Spitzer H 2 S(0) observation was such that the 
slit length (22.3") is along the major axis and the 
slit width (11.1") is perpendicular to the major 
axis. For the short-high IRS Spitzer observations 
of the H 2 S(l) and S(2) lines, the slit is rotated, 
with the slit length (11.3") oriented perpendicular 
to the major axis of NGC 891 and the slit width 



(4.7") along the plane of the galaxy. Thus perpen- 
dicular to the plane of the galaxy the slits cover 
about the same area. This width is also nearly 
identical to the FWHM of the ZEUS beam. How- 
ever, there is a gap of 30" between the individual 
positions of the Spitzer observations along the ma- 
jor axis. Hence, to convert the Spitzer H 2 obser- 
vations to the ZEUS beam size we have to assume 
a distribution of the H 2 emission. 

In our analysis, we also include the CO(3-2) 
line, which is available in the literature (Mauers- 



berger et al.|1999||Dumke et al.|2001[ |Bayet et aT 



2006 ) . Unfortunately, the line intensities and flux 



densities are not consistent. Despite very similar 
beam sizes, the reported values vary by a factor 
of about 6.5. For our analysis we use the inten- 
sity given by Mauersberger et al. ( 1999[ ) ( Table |2j, 
which are about a factor of three less than the 
largest reported intensity. 



The CO(l-O) (e.g. Scoville et al. 1993 Sofue 



fc Na kai 1993 | and CO( 3-2) ( |Mauersberger eTai 



1999||Dumke et al.|200l"| ) emission shows a plateau 
of extended emission with strong central peak and 
local maxima roughly 1' to the North-East and 
South- West. The H 2 Spitzer observations also 
show extended emission with a central peak, a 
trough at the next sampling positions 30" away 
from the center, local maxima at the second sam- 
pling position at 60" from the center on either 
side, and a decrease in the following sampling po- 
sitions. Thus, within the limited sampling of the 
Spitzer observations the distribution of the H 2 and 



CO emission appear similar. Dumke et al. (2001) 



estimated a source size of about 17" x 8" (de- 
convolved) for the central gas concentration from 
their CO(3-2) map. The profile of the CO(l- 
0) emission along the major axis is similar, but 
slightly larger, showing a FHWM of about 25" 



above the plateau (Fig. 2 and 3 in Scoville et al 
1993} |Sofue fc Nakai||1993| , respectively). It is 



plausible that the higher excitation CO emission 
is more concentrated toward the center than the 
CO(1-0) emission. Here, we use the source size 
estimated from the CO(3-2) emission to convert 
the H 2 Spitzer observations and the CO (7-6) and 
[C I] intensities to the 11" beam size of our ZEUS 
CO(6-5) observations, resulting in conversion fac- 
tors of 1.4 for the H 2 S(0) line and 0.9 for the H 2 
S(l) and S(2) lines. The converted intensities are 
given in Table [5] 
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To model the shock parameters we use a dif- 
ferent method than |SCB10| They compared the 
H 2 emission with the far-infrared ISO lines on the 
scale of the large ISO beam, determined the PDR 
model that best fits the data (biased toward the 
ISO far-IR lines), and attributed the remaining 
H2 emission to shock-excited, micro-turbulent gas. 
Since the ZEUS observations cover a much smaller 
region we can not use the ISO observations to 
constrain the line emission fraction. Scaling the 
ISO observations to the much smaller beam size 
of ZEUS is subject to large uncertainties without 
knowing the morphology of the gas-phase locked 
in PDRs. While a reasonable scaling factor could 
be inferred for some far-IR lines, it would be ques- 
tionable for others. Several observations with the 
Kuiper Airborne Observatory (KAO) have found 
that the distribution of the [C II] and CO (1-0) 



we use the observed morphology of the submillime- 
ter continuum. Serabyn et al. (19991 obtained a 



emission are correlated in galaxies (cf. Crawford 
et al.|1985 Stacey et al.|1991 ) . It is therefore rea- 



sonable to use the observed CO(l-O) morphology 
in NGC 891 to scale the [C II] ISO observation to 
the smaller ZEUS beam. However, no such corre- 
lation has been shown for the [O I] lines. In our 
approach we therefore started by only using the 
higher spatial resolution observations, the CO(7- 
6), (6-5), and CO(3-2) transition together with the 
the H 2 S(0), S(l) and S(2) lines, and applied a con- 
strained least squares algorithm to simultaneously 
fit the individual line emission fraction, beam fill- 
ing factor, and the shock parameters for just a C- 
shock and a J-shock model. Unfortunately, trying 
to fit all these parameter to just a shock model 
results in an under-constrained problem and no 
unique solution is possible. This is apparent in un- 
reasonably small x 2 values for a variety of shock 
parameters, line fractions, and beam filling fac- 
tors. In addition, the beam filling factor and the 
line fractions are degenerate. 

To obtain useful model fits it is necessary to 
constrain the fraction of the line emission from 
shocks. This can be achieved by introducing and 
simultaneously fit a second gas component for the 
remaining line emission. An obvious choice for this 
component are PDRs, which have been traced in 
NGC 891 through the far- infrared line emission 
on large scales. In our fitting routine for the PDR 
phase we include our [C I] 370 /im line and the 
infrared continuum. 

To scale the far-IR continuum to our 11" beam 



map of NGC 891 in the 350 fxm continuum us- 
ing SHARC on the CSO. They estimated that 
about 4.8% of the continuum flux density from 
the entire galaxy emerges from within the central 
12". Here we assume that the total IR luminosity 
of NGC 891 of L(8 - 1000/x m) = 1.89 x 10 10 L Q 
(RBGS; Sanders et al.||200'3 ; converted to a dis- 



tance of 9.5 Mpc) scales the same. We convert the 
IR luminosity from the central region into an in- 
tensity and multiply that value by 1.09 to scale to 
our 11" FWHM beam. A large fraction of the IR 
continuum emission in NGC 891 likely arises from 
a diffuse dust component and is probably not asso- 
ciated directly with PDRs. Considering the entire 
galaxy as much as 69% of the IR continuum might 



arise from this component (Popescu et al. 20041 



Given that our observation is focused on the cen- 
ter of NGC 891 the percentage might be lower. 
Here we follow [SCBlol and assume a value of 50% 
as the IR- fraction arising from cirrus in NGC 891, 
and allow an uncertainty of 50%. 

We exclude the [C I] 609 ^m observation (Ta- 
ble [2| since the intensity appears unreasonably 
small compared to the [C I] 370 /jm line intensity. 
The ratio of the (converted to 11") [C I] 370/[C I] 
609 intensity is 17.3. In a PDR this ratio increases 
with density and FUV field, and for a density of 
71 = 10 6 cm" 3 and a FUV field of G D = 10 6 (in 
units of the Habing field: 1.6 x 10 3 erg/s/cm 2 ) 



this line ratio is only 7.8 (Kaufman et al. 1999). 
For a PDR density of n = 10 3 cm" 3 and a FUV 
field of G — 10 2 the expected line ratio would 
be 2.4. The CO(3-2) intensity quoted by |Bayet 



et al. (2006) is also the smallest value compared 



with the values reported by Mauersberger et al. 
(|1999|) and |Dumke et al.| ( |2001[ ). It could be pos- 



sible that our ZEUS [C I] 370 ^m intensity is too 
high. However, we have measured the [C I] 370 /xm 
line simultaneously with the CO (7-6) line. If the 
[C I] 370 /xm intensity is too high would also mean 
the CO (7-6) is too high. Both the [C I] lines and 
the cirrus-subtracted IR continuum are expected 
to originate in PDRs. The [C I] 370/IR and [C I] 
609/IR intensity ratios is 2.7x 10~ 4 and 1.6 x 10" 5 , 
respectively. For a PDR density of n — 10 3 cm -3 
and a FUV field of G = 10 2 the expected ra- 
tios are 3.2 x 10~ 4 and 1.3 x 10~ 4 for the [C I] 
370/IR and [C I] 609/IR intensity ratios, respec- 
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tively. Although we can not rule out errors in 
estimating the IR continuum arising from PDRs 
in the central region of NGC 891 or a systematic 
error in our calibration, it appears that the ZEUS 
[C I] 370 /im intensity is consistent with other 
measurement and PDR predictions. A possibility 
for the smaller than expected [C I] 609 /im inten- 
sity is self-absorption by cold foreground gas. The 
[C I] 609 fim transition connects to the ground, 
has an upper level energy equivalent to 24 K, and 
atomic carbon is the fourth most abundant ele- 
ment. Thus, in an edge-on galaxy self-absorption 
of fine-structure line transitions between the low- 
est energy levels of abundant elements is a plau- 
sible scenario. For example, [O I] 63 /jm line 
also connects to the ground, and self-absorption 
by colder foreground gas has also been invoked to 
explain the relatively weak [O I] 63 /im emission 
from galaxies (e.g. Vasta et al.|2010 and references 
therein) . 

We compare the observed to the predicted line 
intensities from the combination of a PDR and 
shock model and apply a least squares algorithm 
to fit the individual line fractions and the beam 
filling factors. The model grid for shocks (Flower 



& Pineau Des Forets 2010) is rather coarse, with 



just two values for the pre-shock densities (2 x 10 4 
cm -3 and 2 x 10 5 cm -3 ) and velocity ranges from 
10-30 km/s for J-shocks and 10-40 km/s for C- 
shocks, with a resolution of 10 km/s. The model 
grid for the PDR parameters covers a large range, 
with densities, log(n/cm~ 3 ), and far-UV fields, 
log(G ), between 1-5, with a step size of 0.25. In 
the fitting algorithm, we consider the model in- 
tensities from a single PDR surface. We also take 
into account that the CO emission is usually op- 



tically thick (e.g. Bradford et al. 2005) and that 



50% of the infrared continuum originates from cir- 
rus in NGC 891. For an externally FUV-irradiated 
molecular cloud PDRs are created on the outer 
surface. Thus, along the line of sight a beam in- 
tercepts two projected PDR surfaces from a sin- 
gle cloud, one at the front of the cloud and one 
at the back of the cloud. An optically thin line 
would be observed from both PDRs, while an op- 
tically thick line would only be detected from the 
front side. To compare the predicted line intensi- 
ties from a single PDR surface with the observa- 
tions we therefore multiply the predicted intensi- 
ties of the optically thin lines by a factor of 2. The 



beam filling factor derived for PDRs is then with 
respect to single PDR surfaces and not for entire 
clouds. We apply this fitting routine to two com- 
binations: a PDR and C-type shock combination 
and a PDR and a J-type shock combination. Note 
that the beam filling factors can be larger than 
unity if the beam intersects many PDR surfaces 
or shocks along the line of sight, with each con- 
tributing to the observed intensities. We do not 
include the far-infrared ISO lines in the fitting al- 
gorithm because of the uncertain scaling factors. 
Instead, we just list the predicted intensities for 
the [O I], [C II], and [C I] 609 /im lines from the 
shock and PDR models. The results are shown in 
Tables [3j |4j and [5] where the predicted intensi- 
ties from shocks and PDRs are already multiplied 
by the appropriate beam filling factor and for the 
PDRs, the optically thin lines are multiplied by an 
additional factor of 2. 

The fit of combined PDR and shock models 
to the data is slightly better for J-type than for 
C-type shocks. However, most of the resulting 
model parameters are similar and the difference 
in x 2 is n °t enough to be conclusive. In both 
cases the best fit is for a pre-shock gas density 
of n = 2 x 10 4 cm -3 , a PDR density of rtpDR. = 
3.2 x 10 3 cm" 3 , and a FUV field of G Q = 100. 
For the J-shock/PDR combination, the remaining 
parameters that best fit the data are a shock ve- 
locity of 20 km/s, a shock beam filling factor of 
0.87, and a PDR filling factor of 1.34. For the C- 
shock/PDR combination the best shock velocity 
is 10 km/s and the beam filling factors are 0.25 
for the C-shock and 1.28 for the PDR component. 
The fraction of the individual line emissions is also 
very similar for both types of shocks. In both 
cases 100% of the CO(7-6) and 99% of the CO(6-5) 
intensity would arise from shock-excited gas. Of 
the CO (3-2) emission, 47% and 32% would arise 
from shocks in the C-shock/PDR model and J- 
shock/PDR model, respectively. For the H2 emis- 
sion, 94% and 93% of the S(2) intensity and 45% 
and 44% of the S(l) intensity would arise from 
shocks in the C-shock/PDR and J-shock/PDR 
model, respectively. For both model combinations 
only about 3% of the H 2 S(0) intensity would orig- 
inate in shocks. The main difference between the 
results of the model combination is in the beam 
filling factor of the shocks. 

The parameter solution for the shock models is 
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very sensitive to a change in observed intensities, 
and more observations would be very helpful to 
provide better constraints. This is mainly due to 
the coarse grid of the shock model. Marginaliza- 
tion of the probability density distribution shows 
that the PDR solutions are strongly peaked with 
a 1 sigma range less than a step-size in the model 
grid for density and far-UV parameters. While the 
marginalized probability density distribution for 
the shock velocities also strongly peak at the best 
fit values with a 1 sigma range less than a step size, 
running the fitting routine with the CO (3-2) inten- 
sity increased by about 10% changes the densities 
of the C-shock/PDR model to a pre-shock density 
of n = 2 x 10 5 cm~ 3 , a shock filling factor of 0.11, 
and a PDR density of npoR = 5.6 x 10 3 cm -3 . 
In contrast, J-shock/PDR solutions are hardly af- 
fected. The fractions of the line emission origi- 
nating from shocks are quite robust. Given the 
uncertainty in the estimate of the IR continuum 
and the observed CO (3-2) intensity, for which we 
assumed a relative error of 50% in the fitting rou- 
tine, the predicted intensities are not better than a 
factor of 2 in most cases. In fact, a change in shock 
velocity by 10 km/s can especially change the pre- 
dicted [O I] intensities originating from shocks by 
an order of magnitude. Nevertheless, except for 
high-velocity J-shocks the [O I] intensities from 
shocks would always be smaller compared to the 
predicted PDR intensities. 

We note that the observed [C I] 609 fim inten- 
sity is too weak by roughly an order of magnitude 
compared to the J-shock/PDR model result and 
by about a factor of about three compared to the 
C-shock/PDR model result. 

Since there arc no high spatial resolution obser- 
vations available for the far-IR lines we can only 
compare the model predictions among each other. 
Including the beam filling factors, the models pre- 
dict a fraction of < 1% of the [O I] 63 /im and 
146 /jm intensity from C-shocks or J-shocks. For 
our best fit shock/PDR model combinations, the 
entire [O I] emission is predicted to originate from 
PDRs. It is nevertheless interesting to notice that 
the predicted intensity of the [O I] 63 /im line 
is similar to the intensity in the H2 S(0) line in 
shocks. However, both intensities are small com- 
pared to the intensity in the higher H2 transitions. 

We now compare the predicted far-IR line emis- 
sion that likely emerges from the central 11" with 



the ISO observations. According to the CO(1-0) 
cut along the major axis of NGC 891 (e.g 



So- 



ftie & Nakai 1993), about 20% of the line emis- 



sion within the 75" (FWHM) central region orig- 
inates from the central 11" (FWHM). If the ISO 
far-IR fine-structure lines scale similarly, as would 
be a reasonable assumption for the [C II] line, then 
the predicted line intensities from our PDR/shock 
model are only about a factor of 2 larger than the 
scaled ISO far-IR lines. 

We can also just fit the PDR model grid to the 
ISO observations of the center of NGC 891. If 
the IR continuum would scale similar to the [C II] 
emission as observed by ISO then about 23% of the 
IR continuum from the entire galaxy would arise 
from the inner 75" (FWHM). As before, we also 
take into account that probably only 50% of the IR 
continuum might arise from PDRs (with the rest 
originating from cirrus) and that the [O I] 63 /im 
line is optically thick. The best PDR fit to the 
observed [C II], and [O I] 63 /im and 146 /mi ISO 
observations and the IR continuum is obtained for 
a gas density of tipdr = 1 X 10 3 cm~ 3 , a FUV field 
of log(G ) = 2.2, and a PDR beam filling factor 
of 0.07. The observed intensities and the inten- 
sities predicted by the PDR model are shown in 
Table [6j Both, the estimated density and far-UV 
field are similar to the predictions we obtained for 
the combined shock/PDR model within the much 
smaller ZEUS beam. In contrast, the beam filling 
factor is much smaller than the one we obtained 
previously for the ZEUS observation, as would be 
expected for the large ISO beam, whose diameter 
is much larger than the thickness of the edge-on 
galactic disc. This small PDR beam filling fac- 
tor within the ISO beam would correspond to a 
filling factor of 3.26 within the ZEUS beam. It 
thus appears that the PDR conditions are simi- 
lar within a factor of a few (< 10) within small 
region covered by ZEUS and the large region cov- 
ered by ISO. If we apply the beam filling factor 
corresponding to the ZEUS beam (3.26) instead 
of the ISO beam filling factor, which translates 
into multiplying the predicted values in Table|6]by 
3.26/0.07, then the predicted intensities from the 
PDR are within factors of 2-3 for the [C I] 370 /xm 
intensity observed with ZEUS and the CO (3-2) 
and H2 S(0) intensities. The PDR model, how- 
ever, fails to predict the mid-J CO intensities by 
more than 2 orders of magnitude and the H2 S(2) 
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emission by about 1 order of magnitude. In fact, 
the large observed CO(7-6)/CO(6-5) intensity ra- 
tio can hardly be modeled by a PDR at all. There 
is no single PDR model that can fit both the CO(7- 
6)/CO(6-5) and CO(7-6)/[C I] 370 /mi intensity 
ratio (Fig. |2| , clearly indicating the presence of an 
additional heating mechanism. 

4. Energy Sources 

We have shown that the warm interstellar 
medium in the nuclear region of NGC 891 is likely 
heated by a combination of FUV radiation and 
shocks, with the mid-J CO transitions probably 
entirely due to shocks. Since NGC 891 is an iso- 
lated galaxy with only a moderate star formation 
rate (SFR) and no active galactic nucleus we do 
not consider shock excitation through powerful 
mechanism that are encountered in ULIRGs or 
galaxy interactions, where strong X-ray emission 
from an active nucleus, nuclear jets, or galaxy- 
wide gas flows likely drive internal or external 
shock fronts that are suggested by their strong H2 
emission. Instead, we have considered moderate 
internal energy sources as heating mechanisms for 
the PDRs, shocks, and turbulence. 

NGC 891 has a modest SFR of 3.8 M /yr 
( Popescu et al.|[2004 ), which could be responsible 
for the PDR excitation of the far-IR and submil- 
limcter fine-structure lines. Sources of shocks and 
turbulence that could excite the H 2 and CO rota- 
tional line emission include cloud-cloud collisions, 
outflows from young stellar objects, or supernovae 
driven outflows. For example, the warm molec- 
ular gas that strongly emits mid-J CO lines in 
the circumnuclear ring in the Milky Way is likely 



high supernova rate (SNRa) of 0.1 yr 1 in the cen- 
tral ~ 100 pc (Bradford et al. 2003). Supernova 



heated by dissipation of MHD turbulence (Brad- 



ford et al. 2005). Cloud-cloud collisions in the 



material that falls into the central region is possi- 
bly the source for the turbulence there. However, 
a 3 pc scale circumnuclear ring, like the one in 
the Milky Way, would not contribute significantly 
within our 500 pc scale ZEUS beam on NGC 891. 

In the center of NGC 253 the mid-J CO emis- 
sion is probably either due to the decay of turbu- 



lent shocks or due to cosmic ray heating (Hailey- 



Dunsheath et al. 2008). Both mechanisms could 



provide enough energy to heat the warm molec- 
ular gas in NGC 253. The enhanced cosmic ray 
flux in the center of NGC 253 is likely due to the 



outflows could also drive turbulence. 



Roussel et al. (2007) have studied H2 emission 



in the SINGS sample, which is comprised mainly 
of "normal" galaxies and some LINERs/Seyferts. 
For the galaxies powered by star formation they 
found a good correlation between the combined 
power in the H 2 S(0), S(l), and S(2) lines and the 
strength of the 7.7 /im PAH emission. In con- 
trast, the H 2 /PAH ratio in the LINER/Seyfert 
galaxies shows a much larger scatter and the ra- 
tio is usually elevated for stronger H 2 emission. 
This suggests that if the combined power in the 
H 2 S(0), S(l), and S(2) lines is above a certain 
threshold, in relationship to the strength of the 
7.7 /im PAH feature, then the H 2 lines are too 
bright to be excited exclusively within PDRs. For 
H 2 sources above this threshold, they suggest al- 
ternative sources of heat ranging from the X-rays 
from a nuclear engine to interactions with nuclear 
outflows via shocks. SCB10 have compared the 



combined H2 rotational lines with the 7.7 /j,m PAH 
feature along the plane of NGC 891 and found 
that the H2/PAH ratios are near the top (thresh- 
old) values for star forming galaxies in the SINGS 
sample for regions within ~ 5 kpc of the nucleus. 
These inner regions contain most (> 80%) of the 
observed H 2 line luminosity for the galaxy so that 
most of the line emission is at the threshold be- 
tween PDR and shock excited origins. However, 
for the outer regions of NGC 891 the H 2 /PAH ra- 
tio appears to rise by a factor of order 2 which 
places these regions more into the shock excited 
regime as defined by Roussel et al. (2007). The 



rise in the H2 /PAH ratio is due to a faster decline 
in the PAH emission compared with the H2 emis- 
sion in the outer galaxy, and might be due to a 
variety of effects including a lessened importance 
of star formation in the outer galaxy, or lowered 
metallicity in the outer regions of NGC 891. This 
later effect is a plausible cause for the enhanced 
H2/PAH ratio seen for the low metallicity dwarf 
galaxy NGC 6822 in the |Roussel et al.| ( |2007| sam- 
ple. 

For our modeled C-type shocks with a veloc- 
ity of 10 km/s and a pre-shock gas density of 
2 x 10 4 cm~ 3 the total expected CO intensity (up 
to J = 20) is 3.2 x 10~ 4 erg/s/cm 2 /sr (Table Al in 



Flower & Pineau Des Forets 2010). The total in 
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tensity of the H2 rotational transitions is a factor 
of 2 higher, and together with CO they provide ba- 
sically the entire radiative cooling for the C-shock. 
Flux energy in the neutral oxygen line is less than 
1% and therefore negligible. About 54% of the 
mechanical flow energy is converted into line ra- 
diation. Applying the derived beam filling factor 
of 0.25 for C-shocks, the total combined luminos- 
ity of H 2 and CO within the 11" beam is then 
2.17 x 10 6 Lq, which is roughly half the mechani- 
cal heating rate. 

For our J-shock solution the total expected CO 
intensity (up to J = 20) is 1.4 x 10~ 3 erg/s/cm 2 /sr 



(Table A2 in Flower fc Pineau Pes Forets 20101, 



a factor of 4.3 higher than for the C-shocks. The 
strongest CO transitions are around J up = 16, so 
the total CO intensity given above is a lower limit. 
The energy in all H2 lines is a factor of about 8 



higher than from all CO lines (Flower & Pineau 



Des Forets 2010). This results in a combined lu- 



minosity of all the H2 and CO lines of 9.7 x 10 7 L 
within the 11" beam and including the shock fill- 
ing factor of 0.87. Both lines account for about 
82% of the cooling, and the combined luminosity 
is more than 40 times larger than for C-shocks. 
H2O would account for an additional 15% cooling, 
and about 97% of the flow energy of this J-shock 
is converted into line radiation. As for the C-type 
shocks cooling via the neutral oxygen lines is neg- 
ligible for J-shocks with the above parameters. 
In NGC 891, the total SNRa is about 0.042 yr" 1 



Strickland et al. 2004). Following Loenen et al 



(120081) and assuming an energy output of lO^ 1 " erg 
per supernova and a transfer efficiency of 10% the 
expected mechanical heating rate is Lsn = 8.34 x 
10 9 x 0.1 x SNRa L Q (eq.4 in |Loenen et al.|2008| ), 
or Lsn = 3.5 x 10 7 L Q . This is a factor of 10 
larger than the required heating rate for the C- 
shocks, but about a factor of 3 less than required 
to explain the J-shocks. 

Another possibility for the heating source are 
outflows from young stellar objects (YSOs). Using 
equation 2 of Loenen et al. (2008), the ejected en- 
ergy rate from YSOs can be estimated as ^|° w ~ 
3 x 10 29 r^ c v^ ow n erg/s, where Ufl ow is the out- 
flow velocity in km/s, n the density of the sur- 
rounding medium in cm~ 3 , and r pc the size of 
the bubble (in parsec) that is affected by the out- 
flow. YSOs only affect the molecular ISM locally, 



et al. 2008) the resulting luminosity is ~ 16 L , 
and thus more than 10 5 YSO's would be required 
within a 11" beam (~ 500 pc) to account for the 
line luminosity in the (C-type) shock-excited gas. 
Given the SFR of 3.8 M /yr flPopescu et al.|2004[ ) 



and a lifetime of the YSO phase of 10 5 — 10 b yr 



(Loenen et al. 2008) the approximate luminosity 
ranges from 6.1 to 61 x 10 6 L , or 2.4 to 24 times 
the line luminosity from C-type shocks. |Lada| 



et al. (2010) suggest a relationship between the 



star-forming rate and the number of YSO based 
on observations of Galactic star-forming clouds: 
N(YSO) = 4 x 10 6 x SFR A^yr" 1 . Applying 
this relationship to the SFR of the entire galaxy 
NGC 891 we would expect about 1.5 x 10 7 YSO. 
NGC 891 is often thought of as an edge-on analog 
to the Milky Way. Withing the central « 400 pc 
of the Milky Way the SFR was estimated to be 



O.O7M yr" 1 (An et al. 2011), which would 



correspond to about 3 x l(r YSO's. If the energy 
transfer efficiency is similar to SN, YSO outflows 
could provide enough energy to drive C-type or 
J-type shock excitation. 

5. Conclusions 

We have observed the CO (6-5), CO (7-6), and 
[C I] 370 /im lines from the nuclear region of 
NGC 891 using our submillimeter grating spec- 
trometer ZEUS on the CSO. To analyze the 
data we have include observations of the CO (3- 



2) ( |Mauersberger et al.||l999| ) and H 2 S(0), S(l), 
S(2) (SCB10), and IR continuum from the litera- 
ture. We find the emission is best explained by a 
combination of PDRs and shocks. While J-shocks 
provide a slightly better fit than C-shocks the type 
of shocks is not conclusive. 

The CO(7-6)/CO(6-5) intensity ratios are large 
for a PDR paradigm and would require extremely 
high far-UV and density values. In addition the 
CO ratio and the CO(7-6)/[C I] 370 /im ratio can 
not be fitted by a single PDR model. However, the 
mid-J CO ratios are consistent with shock-excited 
warm molecular gas, and our ZEUS observations 
together with Spitzer H2 observations can be ex- 
plained by a combined shock/PDR model. This 
result confirms recent studies of the H 2 S(0), S(l), 



so for a cloud with a radius of 0.1 pc (Loenen 



and S(2) transitions with Spitzer (SCB10), which 
required C-shocks to explain the higher H2 rota- 
tional transitions. The best fitting shock models 
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require a a pre-shock gas density of 2 x 10 4 cm -3 
and shock velocities of 10 km/s and 20 km/s and 
shock filling factors of 0.25 and 0.87 for C-shocks 
and J-shocks, respectively. The most likely energy 
sources for the shocks in the center of NGC 891 
are supernovae and outflows from YSOs. While 
both types of shocks fit our data about equally 
well, from energy considerations it is significantly 
easier to provide the requisite mechanical energy 
for C-shock excitation of the lines than it is for 
J-shock excitation. 

We have estimated the fraction of line emission 
originating from shocks and PDRs and find that 
for both types of shocks all of the CO(7-6) and 
99% of the CO (6-5) intensity arises from shock- 
excited gas. For the CO(3-2) emission, about 
47% and 32% would originate from C-shocks or 
J-shocks, respectively. Our model fits also sug- 
gest a larger fraction of the H 2 S(2) (94%) and 
S(l) (45%) emission from shocks than estimated 
by |SCB10| Only about 3% of the H 2 S(0) emission 
would emerge from shock-excited gas. 

The [C I] 370 fim emission arises from PDRs. 
When fitting the combined PDR/shock model to 
the ZEUS and Spitzer data we find a best solution 
for a PDR density of 3.2 x 10 3 cm~ 3 and a FUV 
field of G = 100 for both types of shocks. The 
PDR filling factors are 1.28 and 1.34 for a combi- 
nation with C-shocks and J-shocks, respectively. 
Our best PDR solution is consistent with [O I] 
63 /im and 146 /jm and [C II] observations on the 
scale of the larger ISO beam. Comparing the line 
predictions from our best fitting shock and PDR 
models shows that the [O I] emission arises almost 
entirely from PDRs. 

We thank the CSO staff for their excellent sup- 
port of the ZEUS observations. We also thank 
the anonymous referee whose comments led to a 
great improvement in the data analysis. This work 
was supported by NSF grants AST-0096881, AST- 
0352855, AST-0705256, and AST-0722220, and by 
NASA grants NGT5-50470 and NNG05GK70H. 
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Table 1 
ZEUS Observations 



Line 


Beam 


fT MB dv 


I 


F 




[arcsec] 


[K km/s] 


[erg/s/cm 2 /sr] 


[W/m 2 ] 


CO(6-5) 


11" 


31.2 ±4.0 


1.056E-5 


3.403E-17 


CO(7-6) 


10" 


35.4 ± 7.8 


1.902E-5 


5.066E-17 


[CI] 370/im 


10" 


65.6 ± 7.8 


3.563E-5 


9.490E-17 
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Table 2 

Observations taken from the literature 



Line 


Observation 


Obs.Area 


fT MB dv 




F 






[arcsec] 


[K km/s] 


[erg/s/cm 2 /sr] 


[W/m 2 ] 


H 2 S(0) a 


Spitzcr 


22.3" X 11.1" 




3.39E-5 


1.97E-16 


H 2 S(l) a 


Spitzer 


4.7" x 11.3" 




9.11E-5 


1.14E-16 


H 2 S(2) a 


Spitzer 


4.7" x 11.3" 




4.71E-5 


5.88E-17 


[CI] 609/mi b 


cso 


14.55" 


11.5 ± 1.0 


1.4E-6 


7.9E-18 


CO(3-2) c 


HHT 


21" 


24 ± 2 


1.017E-6 


1.194E-17 



a Extinction corrected H 2 transitions from |SCB10| 
b From Table B.l in |Bayet et"aL] |2006| ) 
c From Table 1 in |Mauersberger et al.| ( |1999} 



12 



Table 3 

Results of the combined C-shock and PDR model fits 





C-Shock 




C-PDR 




0C-sh a 


n (cm -3 ) 


v (km/s) 


</>p D R a n(cm- 3 ) FUV(G ) 


x 2 


0.25 


2.0E+4 


10.0 


1.28 3.16E+3 100 


7.5 



a 0c-sh and (/>pdr are the beam filling factor for the C-shock and the PDR component, respectively. 
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Table 4 

Results of the combined J-shock and PDR model fits 



J-Shock 




J-PDR 




(t>3-sh a n (cm- 3 ) 


v (km/s) 


<W a n(cm- 3 ) FUV(G G ) 


X 2 


0.87 2.0E+4 


20.0 


1.34 3.16E+3 100 


5.0 



a 0j-sh and 0pdr are the beam filling factor for the J-shock and the PDR component, respectively. 
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Table 5 

Observed/converted Intensities within 11" beam, fraction of line emission from shocks 
and predicted intensities of the c-shock/pdr and j-shock/pdr model fits a 



Line 


A>bs 




f(C) b 


C-Shock 


C-PDR 


f(J) b 


J-Shock 


J-PDR 


CO(7-6) 


1.74E- 


5 


1.00 


1.279E-5 


1.508E-8 


1.00 


1.693E-5 


1.580E-8 


CO(6-5) 


1.06E- 


5 


0.99 


1.091E-5 


1.016E-7 


0.99 


9.773E-6 


1.064E-7 


CO(3-2) 


2.24E- 


6 


0.47 


1.469E-6 


1.663E-6 


0.32 


7.774E-7 


1.743E-6 


H 2 S(2) 


4.72E- 


5 


0.94 


3.221E-5 


2.583E-6 


0.93 


3.411E-5 


2.707E-6 


H 2 S(l) 


8.07E- 


5 


0.45 


3.964E-5 


4.987E-5 


0.44 


3.935E-5 


5.226E-5 


H 2 S(0) 


4.17E- 


5 


0.03 


1.338E-6 


2.299E-5 


0.03 


1.224E-6 


2.410E-5 


[OI]63//m 








2.725E-6 


2.809E-4 




1.399E-6 


2.943E-4 


[OI]146/zm 








2.477E-7 


2.788E-5 




3.935E-8 


2.922E-5 


[CI]370Atm 


3.27E- 


5 






3.138E-5 






3.298E-5 


[GI]609/xm 


1.89E- 


6 






9.647E-6 






1.011E-5 


Total IR C 


1.21E- 


1 






6.650E-2 






6.969E-2 


[CII] 










6.106E-4 






6.399E-4 



a To compare model values with observed values we assumed single PDR surfaces in the model fits. The 
model values of optically thin lines are multiplied by a factor of 2, which assumes that optically thin emission 
arises from twice as many PDR surfaces as optically thick emission. Optically thin lines are: all H2 lines, 
[O I] 146/im, both [C I] lines and IR continuum. The predicted model intensities given in the table are also 
multiplied by the corresponding beam filling factor. Intensities are in units: erg/s/cm 2 /sr. 

b Fraction of observed intensity arising from either C-shock or J-shock. Note that the shock and PDR 
intensities given are values from a model grid with limited resolution and thus don't have to mirror the exact 
intensity fractions for the observed intensities as derived by the fitting algorithm. 



C IR (8-1000/im) continuum scaled according to the 350 /im continuum SHARC map|Serabyn et ah (19991 



see text. About 50% of the infrared continuum likely originates from cirrus clouds in NGC 891. The value of 
the observed/converted IR intensity is not corrected for the cirrus emission, yet. However, model intensities 
were fitted to the corrected IR intensity. 
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Table 6 

Observed intensities compared with predicted PDR intensities within ISO beam a 



Line 


Observed 


PDR 


[CII] b 


3.586E-5 


3.639E-5 


[01] 


1.377E-5 


1.437E-5 


[01] 146/zm 


1.802E-6 


1.561E-6 


Total IR C 


6.266E-3 


6.525E-3 


[CI] 370/zm 




1.249E-6 


[CI] 609(im 




5.031E-7 


CO(7-6) 




3.998E-11 


CO(6-5) 




2.684E-10 


CO(3-2) 




2.695E-8 


H 2 S(2) 




8.030E-8 


H 2 S(l) 




1.947E-6 


H 2 S(0) 




5.511E-7 



a To compare model values with observed values we assumed single PDR surfaces in the model fits. The 
model values of optically thin lines are multiplied by a factor of 2, which assumes that optically thin emission 
arises from twice as many PDR surfaces as optically thick emission. Optically thin lines are: all H 2 lines, 
[O I] 146^m, both [C I] lines and IR continuum. The predicted intensities are also multiplied by the beam 
filling factor (0.07). Intensities are in units: erg/s/cm 2 /sr 

b The observed [C II] intensity is multiplied by 0.7, which assumes that only 70% of the line originates 
from PDRs 

c 50% of the infrared continuum likely arises from cirrus clouds in NGC 891, which is corrected in the 
model values. 
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Fig. 1.— ZEUS [C I] and CO(7-6) (top), and 
CO(6-5) (bottom) spectra. 
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Fig. 2. — Plot of parameter range of PDR 
model (FUV field, G a , in units of Habing field 
1.6 x 10 3 erg/s/cm 2 versus density) depicting the 
observed CO(7-6)/CO(6-5) (dashed) and CO(7- 
6)/[C I] 370 /im (dashed-dotted) intensity ratios. 
The uncertainties are indicated by dotted lines. 
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